Spontaneously Broken Neutral Symmetry in an Ecological System 
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Spontaneous symmetry breaking plays a fundamental role in many areas of condensed matter 
and particle physics. A fundamental problem in ecology is the elucidation of the mechanisms 
responsible for biodiversity and stability. Neutral theory, which makes the simplifying assumption 
that all individuals (such as trees in a tropical forest) -regardless of the species they belong to- 
have the same prospect of reproduction, death, etc., yields gross patterns that are in accord with 
empirical data. We explore the possibility of birth and death rates that depend on the population 
density of species while treating the dynamics in a species-symmetric manner. We demonstrate 
that the dynamical evolution can lead to a stationary state characterized simultaneously by both 
biodiversity and spontaneously broken neutral symmetry. 



Neutral models have been proposed to capture the sta- 
tistical structure of tropical forests [l[ . Even though the 
approach is highly debated Q, the neutral hypothesis 
has led to a general and fundamental framework to study 
both the statics 0] and the dynamics Q of ecosystems us- 
ing general tools borrowed from stochastic processes and 
non-equilibrium statistical mechanics. The fundamental 
assumption of neutral theory [l[ is that within a trophic 
level any individual/organism behaves independently of 
the species it belongs to. In other words, the dynamics 
of the system is unaffected by interchanging/permuting 
species labels of individuals. By using this extremely sim- 
plifying hypothesis many empirically measured statisti- 
cal patterns can be well reproduced [3|— L5| . Going one step 
further, a model can be symmetric -but, strictly speak- 
ing, non-neutral-, a generalization of neutrality where the 
dynamics may depend, for instance, on the local or global 
density of individuals in a community, but no change oc- 
curs on the behavior of a population and on its effects 
on the others in the community upon switching two ar- 
bitrary species' labels 0. In this letter, we address the 
following issues: i) Within a generalized neutral frame- 
work -allowingfor intraspecific density-dependent demo- 
graphic rates [6[- are species able to coexist in a stable 
way up to the temporal scale of speciation which eventu- 
ally averts monodominance and extinction? ii) Can this 
generalized neutral symmetry be spontaneously broken 
so that non-neutral behavior of species can emerge from 
an underlying symmetric dynamics? 

In order to illustrate this, we consider a simple stochas- 
tic model, a variant of the (multi-species) voter model 
0,@|, defined as follows: at every vertex of a regular lat- 
tice of linear size L in d dimensions reside a fixed number 
M of individuals belonging to one of S species. At every 
time step, an individual is picked at random and killed, 
and its place is filled by copying one of its neighbors se- 
lected according to a probabilistic rule to be defined in 



detail below. For illustration, let us consider a generic 
system of S = 4 species and global dispersal where the 
neutral symmetry is not broken (see Fig. la). The frac- 
tion of each species' population fluctuates around the 
same average, 1/4, and is statistically indistinguishable 
from the others. Also, at stationarity the four probabili- 
ties, Pi(n), to find the i-th species with population n are 
identical within statistical errors. In this case, the dy- 
namics of the ecosystem is not changed by any permuta- 
tion of species' labels; however, if each species has its own 
specific parameters for birth, death, dispersal etc., the 
dynamics is no longer symmetric. This explicitly broken 
symmetry makes the previous system of S = 4 species 
behave in a completely different way (Fig. lb). For in- 
stance, if a given species and the remaining are identified 
by distinct sets of parameters, the population fraction of 
one species fluctuates around a given average, 2/5 in this 
case, whereas the other ones fluctuate around a different 
average, 1/5. Even the probabilities Pj's have distinct 
behaviors: three of them are identical and the fourth is 
different as shown in the left inset of Fig. lb. Notice that 
the probability to find a species with n individuals, P(n), 
irrespective of the species identity, has a two-peak struc- 
ture in the non-symmetric case. Unlike the symmetric 
case, a non-symmetric model is necessarily characterized 
by a much larger set of parameters which make the ap- 
proach unsuitable for understanding emergent phenom- 
ena (such as biodiversity). However, we will show in the 
present study that it is possible to define a symmetric 
theory from which non-neutral species' behaviors emerge 
naturally on appropriate temporal scales. This enables 
us to describe species-rich ecosystems with a parsimo- 
nious set of parameters which allows species to coexist 
without the overall symmetry characterizing the model. 
The idea that dynamical symmetry among species can 
be broken is not new in population biology. For instance, 
speciation can be interpreted as a form of bifurcation [§] . 
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However, here we introduce a new concept in community 
ecology which is borrowed from the statistical mechanics 
of phase transitions [T3| , i-e. spontaneously broken neu- 
tral symmetry. As shown in Fig. lc, when the symmetry 
of the model is broken spontaneously, species behave as 
in the non-symmetric case on time intervals shorter than 
a characteristic temporal scale, which will be calculated 
later on. On larger time scales, instead, species' iden- 
tities can be swapped and eventually neutral dynamics 
is recovered. These large temporal scales are also com- 
parable to those at which speciation can occur thereby 
sustaining biodiversity. 

We turn now to the mathematical details of our model. 
Let n" > the population at site x of the ev-th species, 
where a = 1, . . . , S, S being the total number of species. 
Thus J2 a =i n x = M holds for all x and the total number 
of individuals in the whole community is N = ML 2 . In 
the following, we shall also use the alternative variable 
Px = n x/M, the fraction/density of individuals of the a- 
th species at site x. Suppose, that at time t, an individual 
belonging to a certain species 7 and at site x is picked 
at random for removal. Then, call (3 the species' label 
of the individual from one of neighboring sites of x, say 
y, selected to replace it. Note that the dynamics keeps 
the total population per site constant at every time step. 
Thus, a generic n" evolves according to 
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The effective transition rate for this process is propor- 
tional to the population of the 7-th species at site x, 
71%, and to the population of the /3-th species in the cho- 



sen neighboring site, m y . Mathematically, this means 
that the probability of colonization is P(nJ x — > n x T ) = 
K^nJnP. If the proportionality constant, K^y, is cho- 
sen independently of the populations of species at x and 
y and independently of the kind of species involved, we 
get a voter like-model 0, 0] with neutral dynamics (the 
standard voter model has M = 1, i.e. only one indi- 
vidual is allowed to live on each site). In this case, 
regardless of the initial conditions, an infinite size sys- 
tem would inexorably evolve towards a mono-dominant 
state, i.e. an absorbing state where only one of the 5* 
species survives. This is a trivial example of sponta- 
neously broken neutral symmetry. In a more realistic 
perspective, however, different competing effects influ- 
ence species interactions favoring or hampering coloniza- 
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H(, such as, for instance, the Janzen-Connell ef- 
fect in tropical forests [3] , stating that the reproduction 
rate of a given species decreases with its local population 
size, or the Allee effect ,a positive density dependence in 
a small density range 12|>ll3(- Altogether, these effects 



may result in an effective, in general non-linear and non- 
monotonic dependence on the population sizes, 
that we encode in the proportionality constant K^y , now 
dependent, in principle, on the population sizes at both 
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FIG. 1: Example of the evolution of a neutral ecological model 
with 4 species with global dispersal (see main text) for: (a) 
neutral symmetry. All the species are indistinguishable and 
fluctuate around the average value 1/4. In the inset (colors 
are the same as in the main picture) we show the probabili- 
ties Pi(n), and the superposition is perfect within statistical 
errors, and (b) non-symmetric dynamics: species 1 has a dif- 
ferent set of birth and death rates with respect to the other 
three species, and fluctuates around an average density of 
2/5, while the others fluctuate around 1/5. The probability 
Pi(n) differs from the others, as shown in the left inset; in 
the inset on the right, the global probability P(n) is shown, 
c) spontaneously broken neutral symmetry. Here the system 
behaves differently depending on the observation window of 
its evolution: for small time scales, the system appear non- 
symmetric, whereas, for longer time scales, the symmetry is 
recovered. Unlike case (b), all the species show a bimodal 
distribution. The probability P(n) in this case superpose vir- 
tually exactly on the probabilities Pi(n). The total population 
is N = 512 individuals for the case a and b, and N = 2048 
individuals for c. 
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position x and y. However, if the dynamics has to be 
neutral/symmetric then K^y- i) cannot depend explic- 
itly on the species' labels 7 and /?; ii) can at best depend 
only on the densities of species (3 and 7. Indeed, because 
the population of every site is fixed, we obtain the con- 
straint J2 a ^p -yPx = 1 — (Px + Px) which is valid for 
every x and plays an important role in the calculations. 
In order to keep the discussion simple, we consider the 
case = K xy (p^), where p^ represents the density of 
species j3 at y replacing one individual of species 7 at x. 

In order to get some insight into the evolution of the 
ecosystem described above, following the standard ap- 
proach for statistical mechanics systems, let us assume 
infinite dispersal or, equivalently, a well mixed system. 
This assumption - referred to as the mean field limit in 
the physics literature - is useful to simplify the treat- 
ment, while yet capturing the qualitative behavior of the 
model in any finite dimension. In this case, the descrip- 
tion is simple since pj = p 1 for all 7 = 1,...,S, the 
average birth rate of a generic species v is proportional 
to p v (t)K(p v (t)) and the time derivative of J2^=i P^ '(*) 
has to vanish. Thus the evolution equation for the aver- 
age density p v (t) can be derived by a standard Kramers- 
Moyal expansion [l6| of the master equation of our sys- 
tem up to second order: 

P v = Np» [(1 - P ")K{p v ) J2 ? K W\ 

1 (2) 

+ [p v \{\ P v w) + J2 p"K(pn]} 2 t 

where £ = £(t) is a Gaussian white noise J-correlated 
in time. Focusing on the deterministic evolution, we set 
from here on = 0. This is equivalent to neglecting 
fluctuations - of 0(1/ N) smaller than the deterministic 
term - in the analytical treatment. The simulations are 
performed by means of Gillespie's algorithm [l5| consid- 
ering directly the full master equation of the system. 

The neutrality /symmetry of the dynamics is reflected 
in the stationary states obtained when -^p v (t) = 0. Note 
that the drift term on the rhs of Eq.© cannot be de- 
rived from a potential function and therefore the sta- 
tionary states cannot be thought of as minima of an 
analytical function. However, regardless of the form of 
K, there are always S + 1 steady states: one neutral- 
symmetric case, p v = 1/S, v = 1,2,..., S, and S 
mono-dominant situations where only one of the p's is 
1 and the remaining ones are 0. By using local stability 
analysis, one can prove that the mono-dominant states 
are stable only when K(l) > K(0), whereas the con- 
dition K'(jj) < guarantees the stability of the sym- 
metric coexistence. If the function K(z) is linear, Eq. 
([2|) has no other stationary stable solutions. However, 
in a more general non-linear case, new stable solutions 
can show up. It is this non-linearity that allows a spon- 
taneous breaking of the neutral symmetry. The sim- 



plest situation of coexistence within a broken-symmetry 
scenario is obtained when a given species has density 
tp > 1/S and all the other species have the same den- 
sity ( = (1 — <p)/(S — 1) < 1/5, which can occur in S 
different ways. These densities correspond to stationary 
solutions of Eq.© if K(ip) = K(() and are also stable 
when K'{C) < and K'(<p) < -K'(()/(S - 1). We now 
discuss three paradigmatic cases. 

A) K = constant. This corresponds to the classic voter 
model 0,111 (see fig. 2a). The deterministic evolution, 
given by eq.©, is trivial because any initial value of the 
population of each species remains invariant across evolu- 
tion. However, the stochastic dynamics leads to a mono- 
dominant state with only one surviving species, a trivial 
case of spontaneously broken neutral symmetry. For a 
finite system size, the time t(N) to reach one of the 5" 
absorbing states, starting from a random initial condi- 
tion, scales as t(N) ~ iV» where £ = 2, as shown in fig. 
3 (purple line) where log -r(iV) versus log N is plotted. 

B) K(z) = a(b — z) with a, b > 0. This is a more 
interesting case (see fig. 2b) in which the colonization 
ability of a given species at some position decreases as 
its population -at the same position- increases (negative 
density-dependence) and becomes zero when it reaches 
the maximum value b. Therefore, abundant species are 
relatively not as effective in colonizing different regions 
compared to those with small populations. The symmet- 
ric state is the stable stationary state of the determin- 
istic evolution whereas the S mono-dominant states are 
unstable. When the full stochastic dynamics is consid- 
ered, the symmetric stationary state is reached, typically, 
after an initial transient (which depends on initial condi- 
tions). Once the stationary state is reached, it lasts for a 
typical time t(N) ~ exp{«;iV}, as shown in fig. 3 (green 
line) and then the system evolves towards one of the S 
mono-dominant states through a gradual extinction of 
species (observe that this exponential behavior is at vari- 
ance with what happens in the K = constant case where 
t(N) ~ N 1 *). The constant k > depends on the specific 
choice of K(z). The exponential behavior can be easily 
understood focusing on the limiting case of S = 2 species, 
where a description in terms of a potential exists: Intro- 
ducing a density-dependence in the Voter Model dynam- 
ical rule generates an effective potential in the equations 
of motion for p u , v = 1, 2, that in the case of linear KJz) 
discussed above has a minimum for p = 1/2 ( |l7H20| ). 
Thus, applying the well-known Arrhenius law and noting 
that the stochastic term is of order 1/JV smaller than the 
deterministic part (see Eq. [2]), we recover the exponen- 
tial behavior for r(N). For time scales much smaller than 
t(N) or for all times in the infinite size limit, N — > 00, an 
active stationary state exists where all species symmet- 
rically coexist. Therefore, negative density-dependence 
strongly enhances species coexistence. 

We have calculated the relative species abundance 
(RSA) in the steady state, i.e. the probability, P(n), 
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to find a species with population n. The population 
n v {t) of the v-th. species is followed for a sufficiently 
long time and the frequency, f"(ri)An, in each inter- 
val (n, n + An) is recorded and the RSA is obtained as 
P(n) — Ylu=i P u { n )l ' S- In the neutral/symmetric case, 
the P v is independent of v and the corresponding RSA 
is equivalent to those in figure la. Note that at variance 
with the K = constant case -where the RSA is not well 
defined as a consequence of the lack of metastable active 
in the case K(z) 
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a(b — z), (see Fig. la 
obtain a mode, as typically found in the RSA of several 
tropical forests 0, y, 0] and other ecosystems 
C) K(z) has the 'S' shape shown in fig. 2c 
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in order to satisfy the stability conditions for a broken 
symmetry scenario given above (this particular shape 
is for convenience, but it is also valid for K(z) of the 
generic cubic form K(z) = az 3 + bz 2 + cz + d with suit- 
ably chosen coefficients; note that a cubic non-linearity 
in the density-dependence is usually called a Nagumo 
term and is employed to describe populations experi- 
encing the Allee effect [lH). Here the broken-symmetry 
coexistence is the stable stationary state of the deter- 
ministic evolution. Turning on the stochastic dynam- 
ics -after an initial transient- the system reaches one 
of the S stationary states of the deterministic dynamics 
with broken symmetry. Again, on a typical time scale 
t(N) ~ cxpjft'Af} there is gradual extinction of species 
till, one gets a mono-dominant situation. Once more, 
the constant k' > depends on the specific choice of 
K(z). When the system is in a broken-symmetry case, 
the species whose density fluctuates around the average 
(p > 1/S interchanges with one of the S — 1 species fluc- 
tuating around the average £ = (1 — <p)/(S — 1) < 1/5 
on time scales T sw u c h{N) ~ cxp{k s N}. Thus in a fi- 
nite system, N < oo and on a time scale 3> T switc h{N) 
the ecosystem looks neutral/symmetric, i.e. species be- 
have like they were interchangeable. However, for time 
scales <C T SW i tc h(N) or for all times within an infinite 
system, N = oo, the neutral symmetry is spontaneously 
broken and the ecosystem looks as if species were not 
all interchangeable. We have calculated the probabil- 
ity, P v (n), that the z^-th species has population n on a 
time scale smaller than T SW i tc h(N) so as to exhibit the 
characteristics of a broken-symmetry state. The results 
are indistinguishable from those of the case where there 
is no neutral symmetry (Fig lb), in which we run the 
model with two different functions K(z) depending on 
species label: for v = 1 we set K(z) = K\{z) = cti — b\z 
with oi = 3 and b\ = 2, while for v = 2,3,4 we set 
K(z) = Kziz) = d2 — b-2Z with a2 = 2.5 and 62 = 1-5. 
The RSA for the spontaneous symmetry breaking case 
calculated for time scales 3> T SW i tc h(N) is displayed in 
the inset of fig. lc where two peaks appear, showing 
that one of the species behaves differently from the oth- 
ers. In a more general pattern of spontaneous symmetry 
breaking, one can have up to S distinct P^'s produc- 



ing a S'-pcak RSA. Multiple peaks would be resolved in 
the RSA depending on the width and separation of the 
peaks: this scenario is consistent with some recent studies 
on several different ecological communities [23| pointing 
out the possibility of a multimodal distribution of P(n) 
in real systems. 

In conclusion, we have shown that a simple non- 
equilibrium microscopic model for a general S'-species 
ecological community driven by a density-dependent but 
otherwise completely neutral/symmetric dynamics -i.e. 
the dynamic rules governing the stochastic microscopic 
process are insensitive to the species' labels- can show a 
rich and stable heterogeneous biodiversity even at very 
long times. The striking fact is that species can behave 
distinctly by spontaneously breaking the neutral symme- 
try. 




FIG. 2: (Color online) a) (Red solid line) K{z) = 1, corre- 
sponding to the standard Voter Model with many species, b) 
(Green dashed line) K(z) = a(b — z): This definition of the 
function K(z) makes the symmetric state stable against per- 
turbations, and the monodominant states unstable, provided 
a > 0. c) (Blue dotted line) K(z) allowing S stable stationary 
states where the neutral symmetry is spontaneously broken 
by one of the S species. 



le+IO 
le+09 

le+08 

~ le+07 



O 1 le+06 
1 

1 00000 

10000 ... 

1000 




100 



Voter like 
Symmetric Coexistence 
Spontaneously broken symmetry 
t(N)=N 

Log(N) 



1000 



FIG. 3: (Color online) Mean time to extinction t(N) for the 
three different definitions of K(z) in Fig. l[2"|l. calculated in 
the mean field approximation and plotted in Log-Log scale 
varying N from N = 100 to TV = 1000. For K = const. 
(red solid line), t(N) ~ N a with a ~ 2 (red dotted line) as 
expected for a Voter-like model, while the two cases of K = 
b — az (green dashed line), where we chose a = 0.04, b — 1.04, 
and K(z) allowing for a spontaneous breaking of the neutral 
symmetry (blue dotted line) show an exponential behavior 
t(N) ~ e kN . In the inset, we show the same plot in a Log- 
Linear scale, to emphasize the exponential growth. 
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